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In this paper we explore the strong rotation limit of the rotating and stratified Boussinesq equations with 
periodic boundary conditions when the stratification is order one ([Rossby nu mber! Ro = e, [Froude number] 



Fr = O(l), as e — > 0). Using the same framework of Embid & Majdal(ll998l) we show that the slow dynamics 



decouples from the fast. Furthermore, we derive equations for the slow dynamics and their conservation laws. 
The horizontal momentum equations reduce to the two-dimensional Navier-Stokes equations. The equation 
for the vertically averaged vertical velocity includes a term due to the vertical average of the buoyancy. 
The buoyancy equation, the only variable to retain its three-dimensionality, is advected by all three two- 
dimensional slow velocity components. The conservation laws for the slow dynamics include those for the 
two-dimensional Navier Stokes equations and a new conserved quantity that describes dynamics between 
the vertical kinetic energy and the buoyancy. The leading order potential enstrophy is slow while the leading 
order total energy retains both fast and slow dynamics. We also perform forced numerical simulations of the 
rotating Boussinesq equations to demonstrate support for three aspects of the theory in the limit Ro — > 0: 
1) we find the formation and persistence of large-scale columnar Taylor-Proudman flows in the presence of 
O(l) Froude number; after a spin-up time 2) the ratio of the slow total energy to the total energy approaches 
a constant and that at the smallest Rossby numbers that constant approaches one; and 3) the ratio of the slow 
potential enstrophy to the total potential enstrophy also approaches a constant and that at the lowest Rossby 
numbers that constant is one. The results of the numerical simulations indicate that even in the presence of 

the low wave number white noise forcing the dynamics exhibit characteristics of the theory. 
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1. Introduction 



Charneyl (ll948) estimated the orders of magnitude 



For planetary scale rotating and stratified fluid dynamics j 
of different terms in the Euler equations by using typical values of large scale atmospheric motion. From 
these arguments, which included approximate hydrostatic and geostrophic balance, he derived reduced sets 
of equations called the quasi-geostrophic equations (QG) that are widely used in idealized studies of oceanic 
and atmospheric dynamics. In addition to finding equations that govern the large scales Charney also points 
out that the QG eq uations 'filtered out' the incons e quential fast wav e s from the large scale motions. 



In the work of 



Embid & Majdal (1996 



1998); 



Majda & Embidl (I1998h they showed that the QG limit 



([Rossby number] Ro — > 0, [Froude number] Fr — > 0, and Fr/Ro = f/N finite), can also be derived 
from the rotating Boussinesq equations with periodic boundary conditions, Eqs. (I2.8b - (l2.9b . Their asymptotic 
analysis relies on the separation of fast and slow time scales and incorporates the fast waves that were filtered 
out in Charney's analysis. The resulting limiting equations are obtained by averaging over the fast time scale 
and accounts for three-waves interactions of fast and slow waves. Moreover, a rigorous justification of this 
approach was given by a d irect application of Schochet's method of cancellation of oscillations for hyperbolic 



equations, Schochet (1994). In these papers, taking Ro — > corresponds to geostrophic balance and taking 
Fr — > corresponds to hydrostatic balance. When both these parameters go to zero the equations for the 
slow dynamics decouple from the fast and are Charney's QG equations. 

In addition to the quasi-geostrophic regime described above we consider the dynamics of two other dynam- 
ical regimes: 1) the strong stratification limit where the physics is dominated by strong stratification but has 
only weak rotational effects and 2) the strong rotation limit where the physics is dominated by fast rotation 
but is only weakly stratified. 

The first limit, the strong stratificatio n limit, is thought to be important in geophysical fluid dynamics, 



see the review by 



Rilev & Lelond (2000). because it describes flows that occur at length scales between the 



large, quasi-geostrophic scales and the s mall scales where energy is dissipated. Fluid dynamical theory for 



this physical regime has been explored by 



Rilev &Lelongl (2000); 



Rilev & deBruvnKo ps (2003); 



Babin et al. 



(1996, 



1997 



1998 



2002); 



Luw 
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Embid & Maida ( 1998). Parametrically this regime is described by Fr — > 0, Ro = 



O(l), f/N—tO. Here it has also been found that the slow dynamics decouples from the fast and that it leads 
to new equations for the slow dynamics that are not the QG equations derived by Charney. 

One way the slow dynamics of the QG limit differs from the slow dynamics of the strong stratification 
limit is in the role of the zero frequency dispersive waves. To explore this we examine the eigenfrequencies 
of the nondimensional linearized rotating Boussinesq equations, Eqs. d2.8b - d2.9b in the absence of dissipative 
effects, 

WW)-/. (doubleX 

Rot r\K\ 

where Iknl = k 2 + I 2 with k and I the horizontal wave numbers, m the vertical wave number, and |k| = 
k 2 +l 2 +m 2 . There are two kinds of eigenfrequencies. The first kind are the slow vortical modes that have zero 
frequency for all k, and contribute to the potential vorticity. The second kind are the dispersive waves that have 
non-zero frequency but make no contribution to potential vorticity. The latter are the familia r inertia-gravity 



waves t hat are filtered fr om Chamey's QG equations. In the strong stratification limit (see 



(1998); 



Babin et al. 



Embid & Majda 



(119971) ') where Fr — > 0, Ro = 0(1), f/N — > 0, the fast eigenfrequencies are, 
'k H | 



Wi? r (k) = ±- 



w(k) = (double). 



(1.2) 



Again there are two kinds of eigenfrequencies, the slow vortical modes and the fast gravity waves. However, 
here the fast waves contribute to the slow dynamics when kn = 0. That is, one of the wave modes, corre- 
sponding to horizontal averages with zero potential vorticity, has a slow component that resonates with the 
P V bearing vo rtical modes. This mani fests itself in the ver tically sheared horizonta l dynamics (VSHF) mo de 



introduced by 



Smith &Waleffe (2002); 



Maida & Grotd (Il997l) 



lEmbid & M aida ( 1998) and investigated by 1 
In this work, we look at the strong rotation (low Rossby) limit where geostrophic balance dominates but 
the flow is only weakly stratified. These dynamics are parametrically described by the limit Ro —> 0, Fr = 

0(1), f/N^oa. 

This physical regime is releva nt in regions of the deep ocean where stratification is weak but rotational ef- 



fects are dominant. For example. 



Van Haren et al. 



(120051) observe values of N = ± .4/ (2.5 < f/N < oo) 



in the deep Mediterranean Sea and argue that the dynamics in those regions are driven by both weak strat- 
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ification and the horizontal components rotation. In fact, they observe nonhydrostatic motions with vertical 

velocities of the same order of magnitude as the horizontal. Another region of the world where strong rotation 

an d weak stratification have been ob served is in the deep Arctic Ocean. Measurements in the Beaufort Gyre 

by.' 



Timmermans et-al. 



(2007 



,2010) show f /N ps 2 above 2600 meters and //iVwoo between the depths of 
2600 and 3600m. One of the reasons these investigators give for studying the deep Arctic is that in their 2002 
pilot study they discovered the dynamics to be unexpectedly active in the deep ocean. Weak stratification in 



the deep ocean at high latitudes has been noted for the North Atlantic and North Pacific in 



Emer y et al. 



(1984) 



where the y compute mean p rofiles of density and Brunt-Va isala frequency; in the deep waters of the Arctic 



Ocean by 



Jones et al. 



Hevwood et al. 



(2002). Furthermore, warm core 



Woodgate et al. 



(2001) 



d 19951) : and in the Southern Ocean by '. 
eddies with depths of 1000 meters or more have been observed in the Arctic by[ 

In the limit of strong rotation and weak stratification (Ro — > 0, Fr = 0(1), f/N — > oo) the fast eigenfre- 
quencies are, 



u Ro (k) = ±i 



;(k) = (double). 



(1.3) 



Again, there are two kinds of frequencies corresponding to fast inertial waves and slow PV modes. Also, in 
this limit the non PV bearing modes make a contribution to the slow dynamics, but this time it occurs when 
m = 0, which corresponds to vertically averaged d ynamics, which we refe r to as Taylor-Proudman dynamics. 



By using the general framework developed in 



Embid & Majdal (119981) we show that in the low Rossby 



number limit the horizontal and vertical dynamics decouple. In the horizontal the slow equations are the two 
dimensional Navier-Stokes equations along with two conservation laws, the horizontal kinetic en ergy and the 



vertic al vorticity. In the case where the flow is not stratified this is consistent with other work 



Chen et al. 



(2005). The vertically averaged vertical momentum equation is an advection-diffusion like equation that 
couples to the buoyancy through its vertical average. The slow equation for the buoyancy is the only quantity 
that remains fully three dimensional and is advected by all 3 components of the slow velocity. The slow 
equations for the vertical momentum and the buoyancy are coupled and give rise to new conservation laws for 
the coupled w — p dynamics. We also show that the slow modes evolve independently of the fast and that the 
total energy is composed of both slow and fast components, while the potential enstrophy is, to leading order 
in the expansion parameter, purely slow. The reduced equations and their conservation laws are supported 
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by numerical simulations using low wave number forcing. These slow equations are not quasi-geostrophic 

because they are nonhydrostatic. 



2. Nonlocal form of the Boussinesq equations 

The Boussinesq equations for flow moving at a constant rotation about the z axis for vertically stratified 
flow is, 

D , 

—v + / z x v + p L pgz + p Vp = j/Av, (2.1) 
-^■p -bw = nAp, (2.2) 
V • v = 0, (2.3) 

where = + v • V is the material derivative, v = (u, v, w) is the Eulerian velocity, p is the pressure and 
the total density p has been decomposed into p = p Q — bz + p, where p is a constant background reference 
value of the density, b is the density gradient in the vertical, and p is the density fluctuation. We assume b > 
for stable stratification. The parameter / is twice the frame rotation rate, g is the acceleration of gravity, v is 
the kinematic viscosity, and n the diffusion coefficient. 

In order to distinguish the physical mechanisms of fast rotation from weak stratification we use the same 
velocity and length scales for all three components of velocity and in all three dimensions. Therefore we 
nondimensionalize using the following characteristic scales; L is the length scale for the three spatial coor- 
dinates x = (x, y, z), U is the velocity scale, and L/U is the advective time scale. The scale for the density 
fluctuation is bU/N, where N — (gb/ Po) 1 ^ 2 is the Brunt-Vaisala frequency. Then we arrive at the following 
nondimensional quantities, 

Ro=^-, Eu = ^, Re = —, Pr=-, Fr=^ r , (2.4) 
jL pu z v n NL 

where Ro is the Rossby number, Fr is the Froude number, Eu is the Euler number, Re is the Reynolds 

number, and Pr is the Prandtl number. Then the non-dimensional Boussinesq equations for rotating and 
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stratified flow are, 

^v + i-zxv + ^Vp + i-pz=i-Av, (2.5) 

*,_^ w= _!_Ap with V.v = 0. (2.6) 

It is clear that v and p are the evolution variables in Eqs. (12.5t - (12.61 l and that the role of the pressure 
gradient term in the momentum equation is to enforce the incompressibility condition. By eliminating the 
pressure term it is possible to recast the Boussinesq equations exclusively in terms of the evolution variables 
and at the same time to incorporate the incompressibility constraint. This equivalent formulation is however in 
nonlocal form. To write these equations in their nonlocal form take the divergence of the momentum equation 
to find the equation for the pressure, 

^ = VA-(^z.u,-±|-V.(v.Vv)), (2.7) 

where A -1 is the inverse Laplacian operator and u = Vxv=(^,!),w) is the local vorticity. Substitute the 
equation for the pressure into Eqs. (I2.5) - (I2.6) to get the nonlocal form of the Boussinesq equations, 

+ ^z x v + VA-i (ii • u, - - V ■ (v ■ Vv)) + ±fA = i-Av, (2.8) 

15l p -F-r W =Rh-r^ (Z9) 
These equations automatically incorporate the incompressibility condition. Indeed, taking the divergence of 
Eq. ( 12.8b results in (V ■ v) = A( V • v), so that if V ■ v is zero initially, then it remains zero for all time. 

A quantity of fundamental importance is the potential vorticity q = u> a ■ V/5, where u> a = u> + fz. Clearly 
q = q — fb, where the evolution of q = — few + w- Vpis given by Ertel's theorem, 

^ = vAu-V P + K V(Ap)-u; a . (2.10) 

If we scale q with fbFr and ui with fRo then the nondimensional form of q is, 

oz br 

and the nondimensional form of Ertel's equation for q is, 

Dq 1 .dp Ro . Ro . _ Ro _ . 

Dl = ReRr A i Z ~ F~rRe + ~Re ' Vf> + ReTr" ' ^ (2 ' 12) 
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The equations for the global integrated total energy and potential enstrophy are, 



77T" / q z dv = I q^- dv = - 



2 dt 



1 



dt 



RePr 



L 






dz 



dv+0(Ro). 



(2.14) 



The energy equation, Eq. ( 12.131 is independent of the Rossby and Froude number but the enstrophy equation 
Eq. (12.14b depends on the Rossby number, with a leading dissipative term depending on |V^| 2 and the 
remaining contributions involving powers of the Rossby number. 



3. Limiting dynamics for the rapidly rotating Boussinesq equations 

Here we formulate the limiting dynamics for the rapidly rotating Boussinesq equations, i.e. in the limit of 
Ro — > and Fr = 0(1). In doing so we will follow the ap proach developed in in great g enerality by E mbid 
and Majda (1998) and which builds upon earlier work of 



Schochetl (1 1987 



Klainerman & Majda (1981) 



Majda (1984) and 



1994). In fact, the present work complements the work of Embid and Majda which focused 
on the cases where Fr — > with either Ro/ Fr finite or Ro = 0(1). The analysis starts with the recasting of 
the rotating Boussinesq equations in an abstract setting that reveals its key structure. This is followed with the 
asymptotic formulation of the slow dynamics equations in the limit of Ro — > and balanced initial data, i.e. 
without fast inertial waves. Finally, we adapt the theory developed by Embid and Majda to formulate limiting 
dynamics equations in the limit of Ro — > and with fast inertial waves. 

3.1. Abstract framework for the rotating Boussinesq equations 

If we introduce the vector u = (v, p) and let Ro = e, then the rotating Boussinesq equations, Eqs. d2.8t - 
d2.9t become, in abstract operator form, 

du 1 



dt e 
u| t =o = u (x), 



Lpu + Lgu + B(u, u) = Du, 



(3.1) 
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with the operators Lp, Ls, B and D given by, 

u 



LpU 



z a v + VA" 1 ^ 



V 



o 



L s u={Fr)~ 1 
( 



V 



- 10 



B(u,u) 



V 



/ 



Du 



v Vv- VA- x (V-(v Vv)) 
v • Vp 
\ 



\ 



(3.2) 



(3.3) 



(3.4) 



/ 



(3.5) 



(i?e)- x Av 

In the equations above the linear operator L = e~ 1 Lp + Ls splits into a fast piece Lp associated with the 
Rossby number Ro = e and a slow piece Ls associated with the Froude number. It is clear that Eq. (I3.lt 
becomes singular in the limit of e — > 0, and the fast operator Lp will have a dominant role. The remaining 
terms in Eq. ( 13. It are given by the bilinear advective operator B(u, u) and the diffusion operator Du. 

As we mentioned before, if the initial data Uo in Eq. ( 13. U is divergence free, then the solution remains 
divergence free for all time. But in fact more is true, each individual operator Lp, Ls, B and D takes 
solenoidal fields into solenoidal fields. Therefore a natural setting for Eq. (13. U is the Hilbert space X of 
vector fields u = (v, p) in L 2 that are divergence free, V • v = 0, and equipped with the L 2 - norm, 
which is physically equivalent to the total energy, ||u|| 2 = J |v| 2 + p 2 dv. In addition, we assume 2tt- 
periodicity in all the space variables. This choice of boundary conditions considerably simplifies the study of 
Eq. (13. U . particularly the analysis of the operator Lp, and the resulting slow limiting dynamics equations, 
Eq. ( 13.131 >. and the fast wave averaging equations, Eq. ( 13.321 >. The reason for this simplification is the fact 
that the associated eigenfunctions are given explicitly in terms of Fourier modes. In addition, the choice of 
periodic boundary conditions is consistent with the numerical simulations presented in Section 4. Changing 
the domain and the boundary conditions can make the mathematical analysis quite difficult; for example, for 
arbitrary bounded domains it may be impossible to characterize the eigenfunctions of Lp. The choice of an 
infinite domain may change drastically the structure of the null space (slow waves) and the range (fast waves) 
of the operator Lp. 
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One key observation is the fact that the operator Lf (and Ls) is skew-Hermitian in X: for ui and 112 in 



X, 



/ U2L^u 1 G?v=— / (Lp U2)*Ui dv, 
Jv Jv 



(3.6) 



'V Jv 

where u* denotes the conjugate transpose of u. Several important properties follow from this fact. First, Eq. 
( 13. Il l satisfies (in the absence of diffusion) the conservation of energy, Eq. (I2.131 l . This property is shared 
with other important systems in mathematical physics, such as the Euler and the Maxwell equations. Sec- 
ond, according to the Spectral Theorem, skew-Her mitian operators have purely imaginary eigenvalues and an 
orthonormal basis of eigenfunctions, see J 



: |Laxl (120021) . Physically this means that the basic normal mode solu- 
tions of equations represent wave motions. Finally, the null space of Lf, N(Lp), is orthogonal to the range 
of Lp, R(Lp). This can be thought of a consequence of the Spectral Theorem because N(Lp) is spanned 
by the eigenfunctions with zero eigenvalue (slow modes) whereas R(Lp) is spanned by the remaining eigen- 
functions with non-zero eigenvalues (fast modes). This last property will be exploited later in the derivation 
of the slow dynamics equations. 



Next we apply the previous observations to the linear equation 



du 1 T 

— + -L F u = 0, 

at e 



and seek normal mode solutions in the form of harmonic plane waves 



(3.7) 



u(x, t) = r cxp [ik • x — ie 1 w(k) t] 



(3.8) 



where k = (k, I, m) is the wave number, w(k) is the frequency and the purely imaginary number A = iui(k) 
is the eigenvalue of Lp associated with the eigenfunction u = r exp [ik ■ x]. The four eigenfrequencies w(k) 
are given by the dispersion relations 



w(k) = ± m/\k\, oj(k) = (double). 



(3.9) 



Therefore the equations admit slow modes moving on time scales 0(1) when w(k) = and fast waves 
moving on time scales 0(1/ e) when oj(k) ^ 0. The fast waves in this limit are gyroscopic or inertial waves. 
Th ey are waves who owe their existence to the presence of the C oriolis force and were orig inall y described 



by 



Kelvinl (Il880h . Descriptions of these waves can be found in 



LeBlond & Mysakl i 19781) and 



Greenspan 
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1990). Of course, if m = in Eq. ( 13.9b then we only have slow gyroscopic waves. Explicit formulas for the 



eigenvectors r associated with the fast and slow normal modes are given in the appendix. 

3.2. Slow limiting dynamics as Ro — > 

Here we consider the limiting dynamics equations as Ro — > under the assumption that the solution u e (x, t) 
of Eq. (13. Il l evolves only on the slow (advective) time scale. The formal derivation in the context of the 
abstract operator equation, Eq. (I3.lt . is straightforward . We start by assuming that u e (x, t) has the asymptotic 
expansion 

u £ (x,t) =u (x,i)+eu 1 (x,t) + O(e 2 ), (3.10) 
as e — > 0. Plugging u e into Eq. ( 13.11 ) and collecting the contribution of order 0(e _1 ) yields 

L F u° = 0, (3.11) 

that is, u° is in N(Lp ) for all time, or equivalently, u° (x, t) is represented exclusively in terms of slow modes. 
In particular, the initial data Uo(x) = u°(x, 0) is in N(Lp) to leading order in e. The next contribution of 
order O(e ) yields 

^ + Lpu 1 + L s u° + B(u°, u°) - Du° = 0. (3.12) 

The slow limiting dynamics equation is now obtained by projecting Eq. (13.121 1 onto N(Lp) as follows. First 
apply the orthogonal projection P of X onto N(Lp) to both sides of Eq. ( 13.121 . Since u° is in N(Lp) for 
all time so is <9u° /dt, hence P(du° /dt) = du° jdt. In addition, since Lpu 1 is in R(Lp), and N(Lp) is 
orthogonal to R(Lp), then P(Lj?u 1 ) = and any contribution from u 1 is eliminated under the projection. 
Finally, we eliminate the superscript in u° and obtain the slow limiting dynamics equations 

-£+P(L s u + B(u,u)-Du)=Q, (3.13) 
at 

u|t =0 = u (x) e N(L F ). 

Next we notice that, to leading order in e, it is enough for the initial data Uo(x) to be in N{Lp ) to automati- 
cally guarantee that the solution u(x, t) of Eq. ( 13.131 remains in N(Lp) for all time. Indeed, if we integrate 
in time Eq. (13.13b and use the fact that Uo G N(Lp), we conclude that u(t) G N(Lp) for all time. Moreover, 
we will show shortly that in the context of the rotating Boussinesq equations the null space N(Lp) consists 
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precisely of the Taylor-Proudman columnar flows, 



?| dynamics 



Tavlor(1921 



11 



. Therefore we can say that to leading order 



in e, if the initial data is a Taylor-Proudman column (i.e. free of fast inertial waves), then the solution u re- 
mains a Taylor-Proudman column for all time and its evolution is described by the slow dynamics equations. 
When the initial data is a Taylor-Proudman column to leading order in e we say that the flow is in approximate 
Taylor-Proudman balance. Finally we remark that all these formal considerations can be established with full 
mathematical rigor t hrough a direct appl i cation of the general theory of singular limits o f hyperbolic sy stems 



first developed by 



Klainerman & Majdal dl98ll) 



Majda ( 1984), with later additions by 



Schochet(1987) 



To obtain the concrete formulation of the slow dynamics for the Boussinesq equations, Eqs. ( 13.1b - ( 13.5b , 
we need to determine explicitly the null space N(Lp) and its orthogonal projection P. For this purpose it 
is convenient to split vectors and operators into their horizontal and vertical components. Thus, the velocity 
v = (v#,u>) with wh = {u,v), the gradient V = (V#, with Vy = ^-), and the Laplacian 
A = A« + £ with A* = £ + 



The null space N(Lp) of the fast operator Lp in Eq. ( 13.21 ) is characterized by, 

Ox 
d 

m+-A^w = 0, 
oy 

oz 



(3.14) 
(3.15) 
(3.16) 



where u — — is the vertical component of the vorticity and v is incompressible, V • v = 0. From 
Eq. (13.16b it follows that A _1 cj = ip is z-independent of z. Introducing tp back into Eqs. d3.14b - d3.15l i. 
shows that ip is the streamfunction for vj, v# = (— |^), and that v H is incompressible, V# • v# = 
~5x ( — §y ) + ~§y = 0- Since v is incompressible by assumption, then it follows that 757 = 0, i.e. w is also 
z-independent. This shows that N(Lp) consists of Taylor-Proudman column flows, i.e. states u = (v, p) 
with v z-independent and v# incompressible. That no restrictions are imposed upon p is not surprising since 
the fast operator L p only includes those contributions associated with the Rossby number. 
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The orthogonal projection P onto the null space N(Lf) is given by, 



Pu = 



(Vff) a - VffA^tVif • (Vfl), 
Ma 



(3.17) 



v p y 

where (f) z = J 27r f(x, y, z) dz denotes the average in the vertical direction. Therefore the concrete form 
of Eq. ( 13.131 ). the slow limiting dynamics equations for the rotating Boussinesq equations is, 

Ovff 1 

-ST" + v h ■ V ff v ff + V H p = -5- a hv h (3.18) 
at ite 

V ff -v ff =0 (3.19) 

dw I I 

-^r + v ff ' Vh-hj = -—Ahw — (3.20) 
at Re b r 

d i +V - Vp -F-r W -Rk^ (3 ' 21) 
where v = v(x, y, t), p = y, z, t), and V^p = Vff A^ 1 [V» • (vij • VffVjj)]. For brevity of notation 
we omit distinguishing between projected and unprotected variables and instead state that all the variables in 
Eqs. ( 13.1 8b -( f3T2Tb are the result of applying the projection operator, Eq (13.17l >. 

In the slow dynamics the horizontal component of the velocity v# is governed by the 2D Navier-Stokes 
equation. Moreover, v# evolves independently of the vertical velocity w and the density p but it influences the 
dynamics of these variables through the advection terms in Eqs. (13.20b - (13. 21b . The dynamics of the vertical 
velocity w and the density p are strongly coupled. Interestingly, the vertical velocity w evolves according 
to a 2D forced advection-diffusion equation, Eq. (13.20b . with buoyancy forcing given by (p) z , the density 
average in the vertical direction. On the other hand, the evolution of the density p is given by the 3D forced 
advection-diffusion equation, Eq. ( 13.21b . and remains the same as Eq. ( 12.6b in the Boussinesq approximation. 

A consequence of this decoupling of the horizontal velocity from the vertical velocity and the density in 
Eqs. (13.18b - (13.21b is the appearance of additional globally integrated conservation laws which are not present 
in the original Boussinesq equations. First, the 2D Navier-Stokes equation for v# yields the conservation of 
horizontal kinetic energy and enstrophy, 

~4 f W H \ 2 da=-±- [ \V H w H \ 2 da, (3.22) 



2 dt J \ Re j ^ 

[ u 2 da= - — 
2 dt J A Re 



' " f lo 2 da = / |Vrf da, (3.23) 
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where integration is over the horizontal period square A = [0, 27r] 2 . Next, taking the vertical average of the 

equation for p, Eq. (13-2 li t, results in an evolution equation for (p) z , 

§- t (p) z + v ff ■ V H (p) z - ±w = ^Ak(^), (3.24) 

Combining Eq. (13.201 1 and Eq. (13.241 1 results in horizontal conservation laws for the vertical kinetic energy 
and the average vertical potential energy, 



1 d 
2dt 
ld_ 

2dt 



\ w 2 da = -— / w(p) z da-— \V H w\ 2 da, (3.25) 
J A Fr J A Re J A 

f (p) 2 z da=jr[ vo{p) z da - j \V H (p) z \ 2 da, (3.26) 

J A. J A. J A- 



and adding these two equations results in the conservation of vertical energy, 

Additionally, if we define the density fluctuation p = p — (p) z , then from the density equations, Eqs. d3.21| > 
- ( |3.24t , it follows the conservation of potential energy, 

where the volume integral is over the period cube V = [0, 2ir) 3 . Of course, Eq. ( I2.13l l for the conservation 
of total energy still holds for the slow dynamics equations, Eqs. d3.18t - ( 13.211 ). Finally, define the potential 
vorticity q for the slow dynamics equations as the the leading term of the expansion of q given by Eq. ( 12.111 ) 
in powers of Ro = e > 9 = ff = ff ■ Then the equation for the conservation of potential enstrophy is 

The conservation laws above were obtained directly from the slow limiting dynamics equations. In the 
different limiting regimes of str ong stratification and Bur ger number Bu = (Ro/Fr) 2 = 0(1), or strong 



stratification and weak rotation, 



Babin et al. 



(11997 



19981) showed that conserved quantities for the limiting 



equations, like the horizontally averaged buoyancy, correspond to adiabatic invariants for the full Boussinesq 
equations. 

Finally, we remark that although the abstract derivation of the slow limiting dynamics equations, Eq. ( 13.13b , 
is completely general, the concrete form of the equations depend on the explicit calculation of N(Lp), R(Lp) 
and the projection operator P of X onto N(Lp). This calculation is very much dependent on the choice of 
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the domain and the boundary conditions. For example, if we assume the horizontal variables infinite in extent 

and periodicity in the vertical variable, i.e. V = R 2 x [0, 2ir], then N(Lp) is still given by Taylor-Proudman 

columns, Eqs. ( 13.14b - ( 13.161 ). the projection operator by Eq. ( 13.171 ). and the slow limiting dynamics equations 

by Eqs. ( 13.181 ) - ( 13.211 1. On the other hand, if we also assume infinite extent in the vertical variable, i.e. 

V = M 3 , then N(Lp) only admits flows with zero velocity (this is reasonable because the only velocity field 

v that is z-independent and has finite kinetic energy is v = 0) and in this case the slow limiting dynamics 

becomes trivial. However, it is arguable whether a fluid of infinite depth constitutes a reasonable assumption 

in the present context. 

3.3. Fast limiting dynamics as Ro — > 

The derivation of the slow limiting dynamics equations presented in the previous section was based on the 
assumption that the solution evolved only on the slow advective time scale. This assumption is warranted if, 
to leading order in e, the initial data does not include fast inertial waves components. If on the other hand the 
initial data contains inertial waves, then the limiting dynamics equations as Ro —> must be modified to take 
into account the fast inertial waves. The derivation of the fast limiting dynamics for small Rossby number and 
finite Froude number is readily obtained by invoking the very general approach developed in the fundamental 
work of Embid and Majda (1998) on the limiting dynamics of the Boussinesq equations with small Froude 
number and either finite or small Rossby number. In fact, we only have to switch the roles of the fast and 
slow operators Lp and L$ and apply their theory in straightforward fashion. For this reason we are content 
with a summar y of the main points o f their theory and their most relevant conclusions for the present work. 



Following Embid & Majdal (119981) we assume that the solution u e (x, t) of Eq. ( 13. \\ depends on two sep- 
arate time scales, the slow advective time scale t and the fast time scale t = t/e associated with the inertial 
waves. In addition, we assume that for e -C 1 the solution has the asymptotic expansion, 

u e (x, t) = u°(x, t, r)| T=t/e + eu^x, t, r)| T=t/e + 0(e 2 ), (3.30) 

and it is also assumed that u 1 (x, t, r) = o(r), uniformly on < r < T/e, to guarantee the asymptotic 
validity of the expansion. The analysis of Embid and Majda then shows that to leading order in e the solution 
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u c (x, t) of Eq. (13. U is given by, 



u e (x, t) = u°(x, t, r)| T=t/£ + o(l) = e-'/^utx, t) + o(l), 



15 



(3.31) 



where u(x, t) solves a reduced equation obtained by averaging over the fast time variable r, 

lim - / e TLF \L s (e- TLF u) + B(e- rLF u,e- TLF u) ~ D(e- TLF u)] dr = 0, (3.32) 

(Jt T->oo T J 

u(x,t)| t=0 = u (x). 



The fast wave averaging equation, Eq. (13. 32b . supersedes the slow dynamics equation, Eq. (|3T 



inertial waves are present. We also remark that the asymptotic analysis of 



Embid & Majda 



3), w henever 



1998) can be 



justified with compl ete mathematic al rigor via the technique of cancellation of oscillations developed in the 



important paper of lSchochetl (119941) . 

In practice it may be difficult to evaluate the limit over the fast variable r in Eq. ( 13.32t . However, this 
calculation can be performed in the case of periodic boundary conditions. In this case the fast operator Lp 
has an orthonormal basis of periodic eigenfunctions of the form 



u£(x) = e 4k - x r£, 



(3.33) 



where k = (k, I, m) is the wave number, a indicates whether the mode is slow (a = 0) or fast (a — ±1), and 
r£ = (v£, p^). The associated purely imaginary eigenvalue is A£ = iuj£, where the frequency w£ is given 
by Eq. (13.91 . namely co^ 1 = ±m/|k| and w£ = 0. The explicit form of these eigenfunctions is given in the 

Appendix. Next, we expand u(x, t) in terms of the eigenfunctions of Lp, 

l 

= E E "k(*)e ik -X, (3.34) 

k a=-l 

and introduce this expansion into the fast wave averaging equation, Eq. (13.321 1. In order to evaluate the fast 



time averaging in Eq. (13.321 i we observe that 



e Tif (e lk 'X) = e iT <e ikx rg, 



(3.35) 



and also that 



lim — 



1 if u = 



(3.36) 

Oifw ^ 

With these observations we can evaluate all the terms in the limiting fast dynamics equation, Eq. ( 13.32b . and 
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conclude that the Fourier amplitudes erg (i) satisfy the system of differential equations 

+ E *$$$<4<£> + E = E 4 Q '' Q) ^'- (3.37) 

The first sum in Eq. (13.37b comes from averaging the nonlinear advection term B(u, u) with summation over 
the set of three-wave resonant interactions i?£ = {(k', k", a', a")|k' + k" = k, + u;£„ = wg}. The 
second and third sums come from averaging the slow operator Ls and the diffusion operator D respectively, 
with summation over the set — {a'\u>£ — Formulas for the interaction coefficients ^„ 
' a ' and D^ 1 are included in the Appendix. 

The fast dynamics equations for the Fourier amplitudes in Eq. (13 .37b suggest that there is strong interac- 
tion of the fast and slow modes through three-waves interactions, via the quadratic interaction coefficients 
However, it is remarkable that the dynamics of the slow modes (a = 0) proceeds independently 
of the fast modes (a = ±1), making the system of slow and fast modes only weakly coupled. This is because 
all the interaction coefficients B^^,, 1 ^ corresponding to "fast + fast — > slow" interaction are always zero. 
The verification of this fact is given in the Appendix. Of course, the equation for the dynamics of the slow 
modes in Eq. ( 13.37b is nothing more than previously derived slow limiting dynamics equations, Eq. ( 13.13b , 
recast in terms of the Fourier modes. Nevertheless, the slow modes influence the dynamics of the fast modes 
in Eq. ( 13.37b through "fast + slow — > fast " interactions. These conclusions mirror those previously derived 
by Embid and Majda (1998) for the case of small Froude number. 

A remarkable consequence of this weak coupling is that in the absence of dissipation there is conservation 
of energy for the slow and the fast modes separately. The reasoning, which was first provided by Embid and 
Majda (1998) for the case of small Froude number, is reproduced below. First we observe that the solution 
u(x, t) of the fast limiting dynamics, Eq. (13.32b . has a unique orthogonal decomposition in terms of slow and 
a fast component, 

u(x,t) = u s (x,f) + u F (x,t). (3.38) 

In fact, u 5 (x, t) is given explicitly by Eq. ( 13.341 ) with a = and u F (x, t) by Eq. ( 13.341 ) with a = ±1. Next 
we substitute u(x, t) into Eq. ( 13.311 ) and make use of Eq. ( 13.35b to conclude that to leading order in e the 
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solution u c (x, t) has the form 

u e (x, t) = u s (x, t) + e- t/eiF u F (x, t) + o(l). (3.39) 

Moreover, since Lp is a skew-Hermitian operator, then e l l tLF is an unitary operator, and since the eigen- 
functions of Lf are also an orthonormal basis for e t l <LhF (see Lax (2002)), we conclude that 

||u|| 2 = ||u s || 2 + ||e t / £iF u F || 2 = ||u s || 2 + ||u F || 2 , (3.40) 

where ||u|| 2 = J |v| 2 + p 2 dvis twice the total energy (kinetic plus potential). In the absence of dissipation 
Eq. ( 12.131 shows that the Boussinesq equations conserve energy, so that ||u|| 2 is constant in time. On the 
other hand, the slow limiting dynamics equations also conserve energy, so that ||u s || 2 is also constant in 
time. Combining this two facts with Eq. ( 13.401 l, we conclude that ||u F || 2 is constant in time, thus proving 
that the energies of u s and u F are constant separately. This important physical property of the fast limiting 
dynamics equations will be exploited later as a diagnostic tool in the numerical simulations. 

Finally, we remark that the Fourier basis in Eq. (13.331 can be used to study the dependence of other related 
physical quantities on the slow and fast modes. For example, the leading term of the potential vorticity q in 
Eq. ( 12.1 11 1 has the eigenfunction expansion 

9=|^EE mp&Z(t) (3.41) 

k a=— 1 

where p£ is the fourth component of the vector r£ in Eq. (13.331) . Inspection of these vectors in Eqs. ( 1A. 2I > 

- ( 1A. 31 l in the Appendix reveals that 1 = and in consequence the q is composed of slow modes even 
with the presence of fast inertial waves in the fast limiting dynamics. This fact can be used as the starting 
point to prove conservation of potential vorticity, in the weak sense and without dissipation, along the same 
lines developed in Embid and Majda (1998). By contrast, the the eigenfunction expansion of the vertical 
component of the vorticity cj is, 

W = £- |H£E ( k < - KKW e lkx , (3-42) 

y k a=-l 

with and v£ being the first and second components of r£ in Eq. ( 13.33l l. Another inspection of Eqs. ( 1A. 2I > 

- (1A. 3b in the Appendix reveals that fcv£ — lu^ = and we conclude that ui is composed exclusively of fast 
modes. 



18 



Beth A. Wingate, Pedro Embid, Miranda Holmes-Cerfon and Mark A. Taylor 



Run Number Ro Fr f N kf ka ey fctotai To 



1 


1.0 


1.0 


7.0827 


7.0827 


3 


3 


1.0 


256 


20. 


2 


0.3 


1.0 


23.6091 


7.0827 


3 


10.0 


1.0 


256 


20. 


3 


0.2 


1.0 


35.4136 


7.0827 


3 


15.0 


1.0 


256 


20. 


4 


0.1 


1.0 


70.8273 


7.0827 


3 


30.0 


1.0 


256 


20. 


5 


0.08 


1.0 


88.5341 


7.0827 


3 


37.5 


1.0 


256 


20. 


6 


0.05 


1.0 


141.6546 


7.0827 


3 


60.0 


1.0 


256 


20. 



7 0.01 1.0 708.2731 7.0827 3 300.0 1.0 256 60. 

TABLE 1 . This figure tabulates the parameters used in the simulations at 256 3 with low wave number forcing. 

4. Numerical simulations 

The goal of this section is to see if key attributes of the Fr w l,Ro — > limiting dynamics can be 
reproduced in numerical simulations that use low wave number white noise forcing. The three aspects we 
examine are: 1) the columnar structure, 2) the time evolution of the ratio of the Er ^o (slow) total energy 
to the total energy, E, 3) the time evolution of the ratio of the Qr ^q (slow) potential enstrophy to the total 
potential enstrophy, Q. 

For all our simulations, detailed in Table ([T), we use the triply periodic, pseudo-spectral LANL/Sandia 
Direct Numerical Simulation (DNS) code that solves Eq (I2.ll - d2.21 l in a rectangular domains (2D [0, l] 2 or 
3D [0, l] 3 ) with a pseudo-spectral method or 4th order finite differences and a RK4 time stepping scheme. 
The code allows for an arbitrary number of passive scalars and arbitrary aspect ratio grids. The code uses 
MPI for its parallelization along with a 3D domain decomposition (allowing for slab decomposition, pencil 
decomposition or cube decomposition). Since its inception it has been designed for performance on massively 
parallel computers and has excellent scalability, which has been demonstrated on up to 18,000 processors 
running problems as large as 4096 3 (64B grid points). All diagnostics and associated I/O are also implemented 
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with fully scalable algorithms. The parallel FFT at the core of the model is one of the fastest available. It 

uses a custom data transpose algorithm which overlaps inter-process communication with on-processor data 

rearrangement allowing the code to rely exclusively on stride 1, on-processor FFTs (for which the code uses 

FFTW). 

For the code's configuration we turn to a paper bv lSmith & Waleffd (12002b . In that work they not only stud- 



ied the generation of large, slow scales in rotating and stratified flow but they also found that for strongly strat- 
ified flows they recovered the sheet-lik e structures described when Vertically Sheared Horizontal (VSH) dy- 



namics dominates which was discussed 



(2003); 



Babin et al. 



Embid&Maid a(1998): 



Rilev&Lelond (2000); 



Rilev & deBruvnKops 



d 1 997b . In the spirit of those simulations we examine some aspects of the low Rossby 



number limit by 



using a similar simulation set-up. The principal difference between the code set-up of 



Smith & Waleffd (120021) and this work is that instead of using high wave number white noise forcing we 
use low wave number white noise forcing which can be understood by considering the Rossby deformation 
radius, Ld, 

f 



T - N T 
L d - jL, 



or kd = — kr. 

N 1 



(4.1) 



where kf is the peak wave number of the forcing and kd is the wave number of the deformation radius. This 
equation shows that the important horizontal length scales described by kd increase with increasing rotation 
rate assuming N is held fixed (see Table (fTJ). Because of limited resolution we choose kf = 3 for all the 
simulations presented in this paper. 

For the sake of completeness we outline some of the details of the ISmith & Waleffd (120021) code con- 



figuration. First, in order to reduce the effects of viscosity in the intermediate range of scales, a hyper- 
viscosity replaces the Laplacian dissipation used in Eq. ( 12. 11 1. The momentum dissipation is replaced by 
(— l) p+1 ^/ l (V 2 ) p v, and the buoyancy by (— 1) p+1 k/ 1 (V 2 ) p / o, where p — 8 for all the simulations presented 
in this paper. The hyperviscosity, Vh, is, 

'E(k m ,ty 



v h = 2.5 



L2-2p 

tVvrt 1 



(4.2) 



Chasnovl (L1994), where k m is the highest available wave number and E(k m ,t) is the kinetic energy in 



the highest available wave number shell. The hyperdiffusivity used for the buoyancy equation is similar. The 
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random forcing spectrum F(k) is Gaussian with a standard deviation s = 1 and energy input rate e/ = 1 

given by, 

F CT p(-.5( fc -MW (43) 

(2tt) 2 s 

In like manner we also use Rossby and Froude numbers based on the energy input rate e / and the peak wave 
number kf of the forcing, 

(e f (2 7 rfc f ) 2 ) 1 /3 ( (2?rfc )2) i/3 
= — ^ and Ro = ^ — -j^- . (4.4) 

The characteristic scales for time, energy, and potential enstrophy are, 

T = (e f (2nk f fy 1/3 , £ = (e f (2nk f )- 1 f 3 , and Q = ( £/ (2 7 rfc / ) 5 ) 2/3 (4.5) 

These scales are used throughout the numerical section to nondimensionalize time, energy, and potential 
enstrophy. The factors of 2ir appear in the above expressions because the code has a domain of [0, l] 3 . 

Each run is spun up from zero and forced throughout by the low wave number white noise forcing described 
by Eq. ( 14.3b . The parameters N and / remain fixed throughout the simulation. For runs where (Ro < 1) 
columnar structures form during the spin up period. The columns are dynamic and retain their basic columnar 
structure in the quantities Tf and u x throughout the simulation (see Fig. (O). 

The evolution of the kinetic energy, potential energy, and potential enstrophy, nondimensionalized using 
the scales in Eq. ( 14.5b . for selected Rossby numbers are show in Fig. (Q~|). Because the flow is forced in the 
momentum equations the kinetic energy increases with time. After a spin-up time, comparison of the different 
runs shows that the smaller the Rossby number the larger the magnitude of the kinetic energy, and the smaller 
the magnitude of the potential energy. The order of magnitude of the potential enstrophy also decreases with 
decreasing Rossby number. 

4. 1 . Columnar Taylor-P roadman flows 

The classical Taylor-Proudman result is that for constant density flow in geostrophic and hydrostatic balance 
the vertical derivatives of the horizontal and vertical velocities are zero, creating columnar flows. Our theory, 
described by the projection operator, Eq. (13.17) , generalizes the Taylor-Proudman theory to the case when 
the density is not constant and the flow is not in hydrostatic balance. We examine the flow characteristics at 
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FIGURE 1 . To characterize the numerical simulations, this figure shows the time evolution of the kinetic energy, potential 
energy, and potential enstrophy for selected Rossby numbers and Ft — 1. Energy, potential enstrophy, and time have 
been nondimensionalized using the scales in Eq. l |4.5t . The kinetic energy grows in time due to the momentum forcing. 
After a spin-up time, comparison of the different runs shows that the smaller the Rossby number the larger the magnitude 
of the kinetic energy, and the smaller the magnitude of the potential energy. The order of magnitude of the potential 
enstrophy also decreases with decreasing Rossby number. For all the simulations except Ro — 1 columns appear that 
span the depth of the fluid. These columns are dynamic and remain columnar for the duration of the simulation. 

different Rossby numbers to see if Taylor-Proudman flows appear. First we define two averages of u, the x 
component of the horizontal velocity, 



u z = / u(x,y, z,t )2n dz, and u x — I u(x,y,z,t )2irdx, (4.6) 

J L J L 

where L — 1 and t Q is any time after the spin up time. Fig. (|2]l and Fig. (01 show contourplots of u z , which we 
used to identify large-scale horizontal structures, and contourplots ofu x , which we use to identify large-scale 
vertical structures. The left panel of each figure shows contours of u z while the right panel shows contours 
of u x . Fig. (0 shows the case when Ro = 1 (not small). The main structures of this flow are consistent with 
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t = 48. 




FIGURE 2. This figure shows vertical (left side) and horizontal (right side) averages of the horizontal component of the 
velocities for Ro = 1 and Fr = 1 simulations of the full Boussinesq equations. Patterns form on length scales consistent 
with the low wave number forcing of kf — 3 but no columns form. 

kf = 3 forcing scale but never become columnar flows. On the right of Fig. ©, for the case of a smaller 
Rossby number, Ro = .2, the contours reveal a low wave number vortical structure. That this is a columnar 
structure becomes evident by examining the lower right panel which shows the contour plot of u x . There is 



Low Rossby limiting dynamics 



23 



FIGURE 3. This figure shows vertical (left column) and horizontal (right column) averages of the horizontal component of 
the velocities for Ro = .2 and Fr = 1 simulations of the full Boussinesq equations. By time t — 10 columnar structures 
are beginning to form. At t = 20 the dominant columnar structures have formed and remain columnar through the rest of 
the simulation. 
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no similar columnar structure for the Ro = 1 run. Therefore, even with low wave number white noise forcing 

columnar Taylor-Proudman flows spontaneously form if the Rossby number is low enough. 

Though much is known about the formation and instability of vortices in rotating and stratified flow, this is 
the first study that we know of where a constant-in-time white noise forcing creates columnar structures that 
remain columnar throughout the length of the simulation. 

Other studies that discuss the forma t ion of columnar s t ructur e s in rotating turbulence (with no buoyancy) 



can be found in 



Davidson et al. 



(2006); 



Staplehurst et al. 



(2008); 



Sreenivasan & Davidson! (12008). These are 



studies about the creation and evolution of vortices from an initial condition, unlike this work which has a 
constant-in-time white noise forcing, and they do not consider stratification, while this work considers only 
weak strat ification. Despi t e this, t he mechanisms th e y foun d for columnar vortex formation are relevant to this 



study. In 



Davidson et al. 



(2006); 



Staplehurst et al. 



(2008) they begin by considering an initial blob of fluid 



in a rotating environment. They find that when the rotation is strong enough linear wave energy propagation 
is biased along the axes of rotation and that when the columnar vortex appears it remains contained in the 
cylinder that circumscribes the initial blob. When the rotation is weak a centrifugal bursting phenomenon 
prevents any columnar vortices from forming. We can see evidence for this in Fig. <(3j as the wave number 
3 structures elongate and finally form columnar structures. There are two classical laboratory studies of the 
format i on of c olu mnar vortices an d both employ a constant-in-time forcing: the original laboratory studies by 



Tavloild 19211) and 



Davies 



dl972l) . Both investigators studied the dynamics associated with the slow, steady, 
horizontal motion of a s olid obs t acle t hrough a fluid rotating abo ut a vertical axis. One main difference 
between the studies is that I 



Taylorl (119211) did not. Both investigators 



: Davie si (119721) included stratification and ' 
find the formation of columnar structures above the moving topographic feature. While there are considerable 
differences in the kinds of columnar vortices that form, when they did form there was no mention of them 
becoming unstable. 

There is also a large body of literature on the stability of columnar vortices in rotating and stratified flow. 
We restrict this discussion to key work relate d to this study (strong ro tation and weak stratification, i.e. 



nonhydrostatic). We first consider the work of 



Potvlitsin & Peltier! d 1 998b in which they use linear stability 



analysis to study the stability of columnar vortices to three-dimensional perturbations. They use two initial 
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distributions of vorticity: Kelvin-Helmholtz-generated vortices in shear and Kida-like vortices in strain. Their 

conclusion is that an isolated anticyclonic vortex column is strongly destabilized by small values of the 

background rotation, while rapid rotation stabilizes both cyclonic and anticyclonic initial conditions. They 

explain this phenomenon using the Taylor-Proudman theorem. They also discuss the details of the stability 

of anticylonic vortices but since all our columnar vor t ices a re cyclonic we will not describe their other results 

here. The next related paper is IPotylitsin & Peltier! Q2003) where they use direct numerical simulations to 



study the evolution of Kelvin-Helmholtz-generated columnar vortices and verify the ir previous results tha t 



strong rotation stabilizes the columnar vortices. The last work we will examine is 



Otheguv et al. 



(2006). 



This work also uses a linear stability analy sis but begins with an initial condition of co-rotating vortices 



When there is no stratification they find, as 



Potvlitsin & Peltier! d 1998b . that the columnar vortex is elliptic 



unstable and that stronger stratification causes a zig-zag instability. This work is not directly applicable to our 
work because it studies the evolution of an initial condition of co-rotating vortices which our simulations do 
not show but it shows that the evolution of columnar vortices depends on the rotation, stratification, initial 
conditions and forcing. 

4.2. Ratio of slow energy to total energy 

Our theory states that in the absence of viscosity the total energy is composed of both fast and slow dynamics 
but that the ratio of the Ro — > (slow) total energy, Er , to the total energy, E, should go to a constant ( cf. 
the discussion following Eq. ( 13.401 i) 

Er 



E 



->• C for v = 0. 



(4.7) 



In this section we examine the time evolution of Er /E of our numerical simulations for varying Rossby 
numbers. The total energy, E, is given by, 



E = J v 1 -(W\"+p 2 )(2nfdv, 



(4.8) 



where V = 1. The Ro —> (slow) component is computed by projecting the full solution vector, (u, v, w, p) 
onto the null space of the fast operator using Eq. ( 13.17) . The total energy from the slow variables is then 
computed by combining the slow horizontal kinetic energy described by Eq. d3.221 > and the total vertical 
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FIGURE 4. The panel on the left shows the total energy, E (solid line) and the slow total energy, Ea (dashed line). The 
total energy and the slow total energy appear to be on parallel trajectories in time. This is explored more fully in the panel 
on the right which shows the time evolution of the ratio Er /E. As the Rossby number decreases a larger fraction of the 
total energy is slow. Furthermore, the low Rossby number runs appear to be approaching a constant at late times. 

energy described in Eq. (13.271 . 

E Ro = J^(\v H \ 2 + w 2 + (p) 2 z )(27r) 2 da, (4.9) 

where A = 1 and, as mentioned above, all the variables are the result of projecting the full solution vector 
onto the null space of the fast operator. The left panel of Figure shows the evolution of both E (solid 
line) and Er (dashed line) with time. At fixed time, both the slow total energy and the total energy have 
larger amplitudes as the Rossby number decreases. The panel also shows that the total energy and the slow 
total energy appear to maintain a similar ratio as time increases. This is explored more fully in the right panel 
where we plot the ratio Er /E. We find that when Ro = 1 the ratio Er /E maintains approximately the 
same percentage of slow to total energy, but that the slow component is a small fraction of the total. As we 
decrease the Rossby number the time evolution of Er /E shows a gradual increase in the percentage of 
energy that is slow relative to the total. For the smallest Rossby numbers we ran, the ratio quickly increases 
to a value where a substantial fraction of the total energy is its slow component and then gradually increases 
toward a constant close to one. While Er /E is not a constant in time, it is nearly so for the smallest Rossby 
numbers. This suggests that even when there is dissipation and forcing in the system there is a rapid initial 
adjustment of Er /E toward a value that indicates a significant fraction of the total energy is slow, then a 
gradual increase toward a constant close to one. 

T T 

To examine the explicit dependence of Er / E on Rossby number we compute its time average, E Ro /E , 
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FIGURE 5. This figure shows the dependence of the time averaged quantities Er / ' E T and Qjj D /Q T on Ro. As the 
Rossby number decreases the simulations show that a larger fraction of the total energy and potential enstrophy is slow. 

where we compute the average using, 

— 1 Tf 

X T = 7f; ^r^TM* (4.10) 



i=T 



where To is a time immediately after spin-up and Tp is the last time available from the simulation. This 
quantity is shown in Fig. (O where we see the trend that as the Rossby number decreases, the ratio of slow 
total energy to total energy is approaching a constant and that the constant is close to one. This implies that 
at the low Rossby numbers most of the total energy in these simulations is slow. 

4.3. Ratio of slow potential enstrophy to total potential enstrophy 

Our theory shows that in the limit of Ro — > and Fr — 0(1), the potential enstrophy contains only slow 
dynamics (see the discussion below Eq. (3.41). Stated another way, the ratio of the slow potential enstrophy, 
Quo, to the total potential enstrophy, Q, goes to one, 



'Ro 



->• 1 



for 



z/ = 0. 



(4.11) 



Here we examine the time evolution of Qb.o/Q for our numerical simulations for a sequence of decreasing 
Rossby numbers. We compute the total potential enstrophy as, 

1 



Q 



2 q dv 



(4.12) 
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where V — (27r) 3 , and q is the nondimensional potential vorticity described by Eq. ( 12.111 ). The Ro — » 

component of potential enstrophy is computed by projecting the full solution vector onto the null space of the 

fast operator, Eq. (13.17l i. However, since the projection operator does not affect p, the Ro —> component of 

the potential enstrophy can be computed by, 

g -=XKI) 2 " v - (4 - i3) 

The panel on the left of Figure ((6) compares the time evolution of Q (solid lines) with the time evolution of 
Q Ro (dashed lines). As the Rossby number decreases, the gap between the solid and dashed lines decreases, 
indicating that the total enstrophy's composition has a larger slow component. This is explored further in 
the right panel where we show the ratio Qr /Q. As the Rossby number tends to smaller values the slow 
component of the potential enstrophy becomes a larger fraction of the total potential enstrophy. This suggests 
that even with low wave number white noise forcing the potential enstrophy of this limit is dominated by its 
slow component. We also note that as the Rossby number decreases the vorticity is expected to have a larger 
fast component (see Section [33) which means the potential enstrophy defined by Eq. ( 14.121 i can be replaced 
with Eq. fiTBI ). 

T T 

Finally, to see the dependence of the ratio Qr /Q on the Rossby number, we plot the quantity Qr /Q 
where the averages are computed using Eq. (14.101 and are shown in Fig. (0. This figure shows that Qr / Q — > 
1 as Ro decreases the ratio is approaching one but more slowly in Rossby number than that of the energy. 
Dimensionally this means the globally integrated potential enstrophy is dominated by the vertical gradient of 
the buoyancy times the Coriolis parameter while both the vorticity times the Brunt-Vaisala frequency and the 
nonlinear terms have lesser influence. 



5. Summary 

We have e xamined the fast rotatio n limit of the rotating and stratified Boussinesq equations using the 



framework of lEmbid & Maidal (11998b - We have shown that to leading order, the dynamics is composed of 
both fast and slow components but that the slow dynamics evolves independently of the fast. We have also 
derived new equations for the slow dynamics. These include the two dimensional Navier-Stokes equations for 



the slow horizontal velocity, a forced advection-diffusion equation for the vertically averaged vertical velocity, 
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FIGURE 6. The left panel shows the time evolution of the potential enstrophy (solid line) and the slow potential enstro- 
phy (dashed line). As the Rossby number decreases the curves that represent the slow potential enstrophy and the total 
potential enstrophy converge toward each other. This is explored more fully in the right panel which shows the time evo- 
lution of the ratio Qr /Q- For the lowest Rossby number, Ro = .01, the magnitude of the ratio starts to decrease after 
about T — 30 indicating that the slow component of the potential enstrophy decreases. However, at that Rossby number 
the simulations are approaching the limit of what 256 3 simulations can resolve. For the other Rossby numbers, as Ro 
decreases this ratio approaches one, showing that in the limit of fast rotation and order one stratification a large fraction 
of their potential enstrophy is its slow component - even with white noise forcing. 

making the slow dynamics non-hydrostatic, and an equation for the buoyancy which is the only quantity to 
retain its three-dimensional character. In the absence of viscosity and diffusivity these new sets of equations 
conserve a horizontal kinetic energy and vertical vorticity, along with a new conserved quantity that describes 
dynamics between the vertical kinetic energy and the buoyancy. The leading order total energy contains both 
fast and slow dynamics, though their ratio is conserved. The potential energy is found, to leading order, to 
contain only slow dynamics. 

Our numerical simulations, which used low wave number white noise forcing in the momentum equations, 
reveal the emergence of Taylor- Proudman flows as the Rossby number decreases. They also support the theory 
that the ratio of the slow to total energy goes to a constant as the Rossby number is decreased and that the 
constant is close to one. In addition to the energy we also examined the ratio of the slow potential enstrophy 
to the total potential enstrophy. The smallest Rossby number we examined, Ro = .01, was at the limit of the 
scales that could be resolved at this resolution. However, the other Rossby number numbers showed the trend 
that the ratio of the slow potential enstrophy to total potential enstrophy also approaches a constant and that 



30 Beth A. Wingate, Pedro Embid, Miranda Holmes-Cerfon and Mark A. Taylor 

constant trends toward one. These numerical simulations indicate that some of the aspects dynamics derived 

in this paper exist even in the presence of white noise forcing and hyperviscosity. 



Appendix 

Analysis of the fast operator 

The fast operator Lp is defined on the Hilbert space X of 27r-periodic square-integrable vector fields u = 
(v, p) that are divergence free, V • v = 0. In the space X the eigenfunctions of Lp are given by Fourier 
modes of the form u k (x) = e ikx r k where k = (k,l,m) is the wave number and r k = (v k ,p k ) is a 
fixed vector. The divergence free condition reduces to the algebraic constraint v k • k = 0. In terms of 
the Fourier eigenmode u k (x) the eigenvalue equation Lpu^ = A k u k reduces to the algebraic eigenvalue 
problem Lp (k)r k = A k r k , where the matrix symbol Lp(k) is given, for k ^ and k = respectively, by 

/ .. ... .. A 



L F (k) 



1 



kl -(l 2 + m 2 ) 

k 2 + m 2 kl 

— Im km 





, mo) 



/ \ 

0-100 



10 





(A. 1) 



with |k| 2 = k 2 +l 2 +m 2 . The algebraic eigenvalue problem has four purely imaginary eigenvalues, A k = iw k , 



with w k 1 = ±|^| corresponding to fast inertial modes and the double eigenvalue w k = corresponding to 
the slow modes. The associated eigenvectors are given as follows. If \kn\ ^ there are three eigenvectors, 



V2\k H \\k\ 



^ -l\k\ + ikrri^ 
k\k\ + Urn 
^\k H \ 2 



V2\k H \\k\ 



^ l\k\ — ikm ^ 
-fc|k| — ilm 
^\kH\ 2 



V 1 / 



(A. 2) 
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The fourth eigenvector does not satisfy the incompressibility constraint Vk • k = 0. If |kff| = but |k| ^ 0, 

then w k 1 = ±1, w k = 0, and there are four eigenvectors, 



(A. 3) 



but the fourth eigenvector, f k , violates the incompressibility constraint. Finally, if |k| = then there are 
two fast modes associated with Wg 1 = ±1 and two slow modes associated with loq = 0, with the four 
eigenvectors in Eq. dA. 31 now satisfying the incompressibility constraint. Notice that the eigenfunctions are 
normalized and satisfy the symmetry condition r£ = r~£, where the bar stands for complex conjugation. For 
this reason we require that the amplitudes er k (t) in Eq. ( 13.341 i satisfy the condition cr k = (7~ k to ensure that 
u(x, t) in Eq. ( 13.34b is real valued. 
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Formulas for the interaction coefficients 



^(a',a",a) j (a* ,a) 



Here we collect the formulas for the interaction coefficients B) k , ' k „ ' k % , L k 



and D 



{a' ,a) 



which appear in 



Fourier formulation of the limiting fast dynamics equations, Eq. ( 13.371 ). The quadratic interaction coefficient 



(k',k",k) ls given by, 



n ru> ic" 



(vg, • k")r£„ + ( V £„ • k')rg, 



(A. 4) 



'(k',k",k) 2 

With this formula we can verify the claim that the fast limiting dynamics equations for the slow modes is 
independent of the fast modes, i.e. that the quadratic interaction coefficients corresponding to "fast + fast — > 
slow" are zero. Because the formula for the quadratic interaction coefficient in Eq. ( IA. 41 ) is invariant under 
the permutation of a' and a", and k' and k", it is sufficient to check that B^^,fK is zero. The three-wave 
resonance equations for this case is 



k' 4- k" — k — 
k +k -k, 



= 0. 



(A. 5) 



There are three cases to consider. First, if |k#| ^ then r k in Eq. dA. 21 i is orthogonal to r k in both Eq. 
(IA. 21 ) and Eq. dA. 31) . and B)^^,fh is zero . Second, if \k H \ = but |k| ^ then r k in Eq. dA. 31 ) coincides 



with r k in the previous case and the quadratic coefficient is again zero. Finally, if |k| = then there r k in 
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Eq. dA. 31 is either r£ or f£. In this case k' = — k" and, by symmetry, rZ^, = r^,. Direct calculation then 

shows that the third component of (v^., ■ (— k'))rj c , + (vL • k')r^, in Eq. dA. 4b is given by, 

(vi, • (-k'))<; + . k ')< = *|k^i 2 ( v i, + ^) • k' 

= 2i|k^| 2 (-Z'|k'|, fc'|k'|, 0) • (k', I', to') = 0, (A. 6) 

and this implies that the dot product with either r£ or f £ for the quadratic interaction coefficient in Eq. ( |A. 4) 
is again zero. This proves that B^^fh is always zero. 

Next, the linear interaction coefficient il a ' a ' in Eq. (13.371 1 is given by, 



L 



(a' ,ck) 



(r£)*L s (k)r£ 



(A. 7) 



where the matrix symbol Ls(k) associated with the slow operator L$ is given, for k ^ and k = 
respectively, by 







km 
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1 




1 
0-10 



(A. 8) 



Direct calculation of these coefficients with the eigenvectors given in Eqs. dA. 2\ - ( 1A. 3t gives the explicit 
values of L^ ±1} = ±i/y/2 and i[ ±1,0) = tvV 2 " whenk = (k,l,0), L { °- 0) = 1, and 4°' 0) = -1 when 
k = 0, and zero otherwise. 

Finally, the diffusion coefficient D^" ' a ' in Eq. (13.371) is given by, 



D 



(a. .a) 



(A. 9) 



where -D(k) is the diagonal matrix 

D(k)=diag(_^_| k |2 ) _^| k | 2i _^ |k | 2) __^_| k 
and direct calculation with the eigenvectors given in Eqs. (1A. 2\ - dA. 31 ) shows that 



(A. 10) 
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(i,i) 



(-i-i) 
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flePr 



|k| 2 , 



(A. 11) 



and zero otherwise. 
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